From Elasticity to Hypoplasticity: Dynamics of Granular Solids 
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Abstract 

"Granular elasticity," useful for calculating static stress distributions in granular media, is gen- 
eralized by including the effects of slowly moving, deformed grains. The result is a hydrodynamic 
theory for granular solids that agrees well with models from soil mechanics. 
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Granular media has different phases that, in dependence of the grain's ratio of deformation 
to kinetic energy, may loosely be referred to as gaseous, liquid and solid. The first phase is 
relatively well understood: Moving fast and being free most of the time, the grains in the 
gaseous phase have much kinetic, but next to none elastic, energy pQ. In the denser liquid 
phase, say in chute flows, there is less kinetic energy, more deformation, and a rich rheology 
that has been scrutinized recently [2]. In granular statics, with the grains deformed but 
stationary, the energy is all elastic. This state is legitimately referred to as solid because 
static shear stresses are sustained. If granular solid is slowly sheared, the predominant part 
of the energy remains elastic. Yet no theory is capable of accounting for both its statics and 
dynamics, and no picture exists that helps to render its physics transparent. 

Two grains in contact are initially very compliant, because so little material is being 
deformed. As this geometric fact should also hold on larger scales, for many grains, diverging 
compliance at diminishing compression is a basic characteristics of granular solids, and the 
reason it is sensible to abandon the approximation of infinitely rigid grains. Starting from 
this observation, a theory termed GE (for "granular elasticity" ) was constructed to account 
for static granular stress distributions. Taking the energy w as a function of Uij, the elastic 
contribution to the total strain field £y, we specify [3] 

w = VA (£§A 2 + Au*) = BVA (§A 2 + ul/0 , (1) 

with A = — uu.i u 2 s = u®jii®j, u®j = — \uu (The notations: a°- = ciy — |a« 5y and a 2 = 
a>ija>ij with any a^- are employed throughout this paper.) The elastic coefficient B, a measure 
of overall rigidity, is a function of the density. Denoting p g as the granular material's bulk 
density, and e = p g /p — l as the void ratio, we take B = B x (2.17 — e) 2 /[1.3736(l + e)], with 
£>o,£ > two material constants. The elastic energy w contributes 7T^ = —dw/duij to the 
total stress a^. And since the elastic stress is the only contribution in statics, force balance 
reads VjCy = V ' = pG{. This was solved for three classical cases: silos, sand piles and 
granular sheets under a point load, resulting in rather satisfactory agreement to experiments, 
see [I]. Moreover, the energy w (with P = |7T«) is convex only for tt s /P < a/2/£, implying 
no elastic solution is stable beyond it. Identifying this as the yield surface gives £ ~ 5/3 for 
natural sand. 

When granular solid is being slowly sheared, we must expect a qualitative change of its 
behavior: In addition to moving with the large-scaled velocity Vi, the grains also move and 
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slip in deviation of it - implying a small but finite granular temperature T g . As a result, some 
of the grains are temporarily unjammed, with enough time to decrease their deformation. 
This depletes the elastic energy and relaxes the static stress. Stress relaxation is typical 
of viscoelastic systems such as polymers. Granular media are similar, but they possess a 
relaxation rate that vanishes with T g . This is the reason they return to perfect elasticity 
when stationary. The basic physics of granular solids, viscoelasticity at finite T g , is in fact 
epitomized by a sand pile, which holds its shape when unperturbed, but fails to do so when 
tapped. A set of differential equations termed granular solid hydrodynamics (gsh) is derived 
consistently below starting from GE, with this simple physics as the only additional input. 
Conservation of density and momentum always holds, 

§- t P + Vi(pvi) = 0, §t(pvi) + Vj(erij + pViVj) = pGi, (2) 

where Gi is the gravitational constant. In granular gas or liquid, the stress cr^ has the same 
structure as in the Navier-Stokes equation, though the viscosity is a function of the shear. 
In granular solid, the stress is not usually taken to be given in a closed form. Instead, 
constitutive relations are employed. These relate the temporal derivatives of stress and 
strain, giving ^cjjj as a function of = |(VjVj + Vjfj) and density (where is often 
replaced by an objective derivative say from Jaumann). 

Hypoplasticity, or hpm (for hypoplastic model), is a modern, well- verified, yet compara- 
tively simple theory of soil mechanics [5]. It is quite realistic in the above specified regime 
of solid dynamics, though less appropriate for determining static stress distributions. The 
starting point is the rate-independent constitutive relation, 

m a v = HijkfVn + Kj\/v ek v ek + e{y u f , (3) 

where the coefficients H^w, Ay, e are functions of &ij,p, specified using experimental data 
mainly from triaxial apparatus. Great efforts are invested in finding accurate expressions 
for them, of which a recent set [3] is e = 1/3, 

H ijk £ = f (F 2 S ik 5 je + a 2 a ij a ke /al n ) , (4) 
Kj = aff d F (aij + (7° ) /a nn , (5) 

where [with a = 2.76, h s = 1600 MPa, = 0.44ej, e c = 0.85ej, e~ l = exp (<t«//i s ) ' 19 , e the 
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If GSH as derived below from the idea given above reduces to hpm under certain con- 
ditions, we would have, on one hand, captured valuable insights into the physics of this 
field-tested theory, understood its range of validity, how to widen it by appropriate modi- 
fications, and on the other hand, obtained a broadside verification of GSH, along with the 
physical picture embedded in it. As we shall see, GSH indeed reduces to Eq (J3J) for a station- 
ary T g , with Hyik, A^, e given in terms of = —d 2 w/duijduki (known from ge) and four 
new scalars [combinations of transport coefficients such as viscosities and stress relaxation 
rates, see Eq (17)]. Although quite different from Eqs (4[5), the new Hijik,Aij,e yield very 
similar accounts in all cases we have considered. 

A large part of GSH may be duplicated from the hydrodynamic theory of transient elastic- 
ity, constructed to describe polymers [Bj. This theory accounts for any system in which both 
the elastic energy and stress relax, irrespective how this happens microscopically - whether 
due to polymer strands disentangling, or the grains unjamming. (A formal and rather more 
detailed derivation of GSH can be found in an accompanying paper [7j.) The stress cry and 
the elastic strain determined by 



a 



it 



fly -Vij, {§i + v k V k ) = Vij + Xij, (6) 



where itij = —dw/diiij is the elastic stress and = \(ViVj + VjV{). a® and are the 
irreversible contributions, given by Onsager relations that connect the "currents," a^,Xij, 
to the "forces," Vij,iiij, 
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L> ~ (V + Vg)Vij + (C + (g)8ijV a + aTTij, (7) 



X^ = -avij + (3^ + Pxdijitu (8) 
= (xv ij - - ^SijUu- (9) 

The coefficients 77, £, r) g , ( g > in 0® are viscosities, see below for their differences. Calculat- 
ing §- t (Jij as in Eq O), they all vanish for steady velocities, = 0. The term X^, accounting 
for the relaxation of the elastic strain u^, is rather more consequential. Eq (|9j is obtained 
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by taking the derivative of Eq Q, vr^ = -dw/dmj = \fK{B/^8 ij -2Au%)+A{u 2 s /2VK)8 ij . 
So the relaxation times are given as 1/r = 2(3 AVA, 1/ri = 3(3i\/~K{B + \Au 2 s / 'A 2 ). The 
coefficient a is a cross coefficient of the Onsager matrix. It is taken as a scalar for simplicity. 

In principle, the transport coefficients rj, i] g , (, ( g , r, r 1; a are functions of the ther- 
modynamic variables: density, temperature and the elastic strain Uy. We shall, again for 
simplicity, assume that they are strain-independent, while noting three points: (1) Constant 
t, T\ implies strain-dependent (3,(3\. Choosing the former as constant and not the latter, 
the trace and traceless part of §~ t Uij are decoupled. (2) As discussed above, 1/r, 1/ti vanish 
with T g . So the obvious and simplest assumption is 

1/r = XT g , 1/n = XxT g , (10) 

with A, Ai, ti/t = A/Ai possibly functions of the density, but independent from stress and 
T g . (3) Being reactive, a is not restricted in its magnitude. It may stay constant while 
1/r, 1/tl vary - though it must eventually vanish for 1/r, \ jr\ — > 0, as a = in statics. 

The above hydrodynamic theory is closed if we amend it with an equation of motion for 
T g . In thermodynamics, the energy change dw from all microscopic, implicit variables is 
subsumed as Tds, with s the entropy and T = dw/ds its conjugate variable. From this, we 
divide out the kinetic energy of granular random motion, executed by the grains in devia- 
tion from the ordered, large-scale motion, and denote it as T g ds g , calling s g , T g = dw/ds g 
granular entropy and temperature. In other words, we consider two heat reservoirs, the first 
containing the energy of granular random motion, the second the rest of all microscopic 
degrees of freedom, especially phonons. In equilibrium, T g = T, and s g is part of s. But 
when the granular system is being tapped or sheared, and T g is many orders of magnitude 
larger than T, then this leaky, intermediary heat reservoir produces physics in its own right. 
Taking s g as the part of the entropy accounting for the granular kinetic energy, our def- 
inition is fairly close to the entropy of granular gas pQ, though its functional dependence 
is probably dominated by the effect of excluded volumes. The entropy s, on the other 
hand, is closer to the so-called "configurational entropy," [8] (see section 6 of the first of jl] 
for a discussion of their relationship). The balance equations are t|s + Vk(svk) = R/T, 
§- t s g + V k {s g v k ) = Rg/Tg, where 

R = V v 2 s + Cv 2 u + (3tt 2 s + (3^1 + 7 T 2 , (11) 
R 9 = Vgvl + Cgvl-lT 2 . (12) 
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The first four terms in the entropy production R are the usual contributions from shear flow 
and stress relaxation, as given by transient elasticity. The first two terms of R g account 
analogously for shear excitation of random motion. The term r fT^ (with 7 > 0) describes 
how the kinetic energy of random motion seeps from s g into s. (Diffusion of T, T g are easily 
included when needed.) 

With Eqs fllj |2j [Hj [7[ [9| [T0| |l2j|llj ), GSH is complete. It especially contains the equilibrium 
case, Oij = 7Ty, in which the dissipative fields vanish, a®,Xij = 0. Off equilibrium, these 

§l<Jij assuming J^i 



two fields are finite, and we calculate JtOVj assuming Jj-Uj = 0, from Eqs (pi [7l [9 



(Tij = (1 - a)§- t Kij = (1 - a)M ijH j^u u 



(1 - a)M ijM [(l - a)v M - \u\f, - ^5 ke u u ]. (13) 

As mentioned above, the energy w looses its convexity at tt s /P = a/6/5, and no static, 
elastic solution is possible beyond this ratio. Therefore, it was identified as yield. Given 



Eq (13), the same identification holds dynamically: The loss of convexity implies that one 
of the six eigenvalues of M^i = —d 2 w / 'du^du^ (written as a 6 x 6 matrix) vanishes at this 
point, and a strain rate along the associated direction yields vanishing stress rate. 

For R g = 0, when s g is being produced and leaking at the same rate, we have a stationary 
T g , given as 

T 9 = ^fikh\Jv 2 s + {t g h g )vl. (14) 



Inserting Eqs ( TofM ) into (13), we retrieve Eq with 

Hijki = (1 - a) 2 M ijki , e = (g/%, (15) 
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[1 - a)M m [(r/Ti)ASkt - u M )X^fnJl- (16) 



hpm has 43 free parameters (36+6+1 for Hijke, Ajj, e), all functions of the stress and density. 
Expressed as here, the stress and density dependence are essentially determined by Mijkt 
that (with £ = 5/3 and Bq = 8500 MPa) is a known quantity |4j. For the four free constants, 
we take 

1 - a = 0.22, — = 0.09, ^ = 0.33, A A j^- = 114, (17) 
n Vg V 7 

to be realistic choices, as these numbers yield satisfactory agreement with hpm. Their 
significance are: CgMg = 0-33 implies shear flows are three times as effective in creating T g 
as compressional flows. t/t\ = 0.09 means, plausibly, that the relaxation rate of shear stress 




FIG. 1: The stress changes d<7i, do"3, calculated using GSH (granular solid hydrodynamics) and hpm 
(hypoplastic model), for given strain rate starting from different points (depicted as crosses) in the 
stress space spanned by a\,a 3 . The strain rate has varying directions but a constant amplitude, 
y/2vf + u|, such that the applied strain changes form circles around each cross (not shown). 

is ten times higher than that of pressure. For a purely elastic system, Eq ^ is replaced by 
jL&ij = MijikVik- Therefore, the factor (1 — a) 2 accounts for an overall, dynamic softening of 
the static compliance tensor My^, a known effect in soil mechanics [9]. Finally, A controls 
the stress relaxation rate for given T g , and \Z%~py how well shear flow excites T g . Together, 
Ay rjg/^ = 114 determines the relative weight of plastic versus reactive response. (Note 
\Aij\f\Hijki\ ~ \u1 e \ ■ 114/(1 — a) is, for around 10~ 3 , of order unity.) 

Next, we compare Eqs (jT5j [l6j) to (|4j |5| in their results with respect to "response en- 
velopes," a standard test in soil mechanics for rating constitutive relations [3]. Axial sym- 
metry of the triaxial geometry is assumed, with &ij,Vij diagonal, and <j\ = a xx = a yy , 
^3 = &zz, vi = v xx = v yy , v 3 = v zz , P = \o x + !<t 3 , q = a 3 - a u a 2 s = |g 2 , d7 = (vi - v 3 )dt, 
de = — (2^1+^3 )dt. Starting from a point in the stress space (spanned by ai,a 3 in Fig 1 and 
a s , P in Fig 2), one deforms the system for a constant time dt, at given strain or stress rates, 
while recording the change in the conjugate quantity. Varying the direction, the input is a 
circle around the starting point, but the response envelopes show deformation characteristic 
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of the system, or the constitutive relation to be rated. Fig 1 and 2 show respectively the 
responding stress and strain envelopes, for the void ratio e = 0.66, calculated using GSH and 
hpm. The similarity in stress-dependence and anisotropy is obvious. 

In Fig 3, one strain envelope is blown up for a more detailed comparison, using the 
extended version of response envelope as given in [TO]. Here, the applied stress rate is 
reversed at halftime, such that the system returns to the starting point in stress space at 
the end. The responding strain change, depicted as deflected, straight dotted lines, does 
not return to the origin. Both GSH and hpm predict that the end points from all angles of 
stress changes (some of the angles are given at the deflection points) form a straight line 
OA. (Instead of a line, a narrow ellipse is reported in the 2D-simulation of [10J. This may 
be a result of the fact that the stationarity of T g is briefly violated when the stress rate is 
reversed, during which the system is rather less plastic.) OA's angle a in strain space is 
usually referred to as the "flow direction," while the direction in stress space, along which 
the plastic deformation is largest (with the strain starting at O and ending at A) is called 




FIG. 2: The change in strain d7, de for given stress rate starting from different points in the stress 
space, spanned by a s , P. The amplitude of the stress rate \J dP 2 + dq 2 is constant. See Fig 3 for 
an explanation of the "flow direction." 
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FIG. 3: A pair of blown-up strain envelopes from Fig 2, with the starting point O at P = 0.2, 
a s = 0.16 MPa. The stress rate is reversed at halftime, and the stress returns to the origin O at 
the end. The strain (depicted as dotted lines) gets deflected, and ends somewhere along OA, a 
straight line for both GSH and HPM. a, the angle of OA, is called the "flow direction;" <p is the 
"yield direction," along which the plastic flow is maximal, with the strain ending at A. 

the "yield direction" 0. Since they are not equal, the flow rule is "non-associated." In Fig 4, 
the flow direction cr, the yield direction <ft, and the maximal plastic strain (the length of 
OA), are displayed as functions of <J S /P, with P = 0.2 MPa. Again, the similarity between 
both theories is obvious. 

We take all this to be a preliminary confirmation for the basic idea of slowly sheared 
granular solids being viscoelastic, and also for GSH as the appropriate hydrodynamic theory. 
Next, it should be interesting to use GSH for circumstances, in which T g is not stationary 
and the stress rate possesses a more complicated form than that given by Eqs 
These include especially sudden changes in the direction of the strain rate j9], such as in 
cyclic loading or sound propagation. Also, one needs to understand whether GSH holds at 
transitions from granular solid to liquid, from t>jj = to Vy ^ for a stationary stress, 
■ir&ij = 0, in phenomena such as shear-banding. 
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P = 0.2 MPa 

GSH 

HPM 




FIG. 4: Yield direction, flow direction, and the maximal plastic strain (length of OA), versus cr s /P, 
for P = 0.2 MPa, calculated employing GSH and hpm, respectively. 
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